Asymptotic behavior of regularized shock solutions in coating flows 



D. Badali^, M. Chugunova^, D.E. Pelinovsky^, and S. Pollack^ 
^ Department of Chemical and Physical Science, University of Toronto at Mississauga, 
Mississauga, Ontario, Canada, L5L 1C6 
^ Department of Mathematics, University of Toronto, Toronto, Ontario, Canada, M5S lAl 
Department of Mathematics, McMaster University, Hamilton, Ontario, Canada, L8S 4K1 
^ Department of Mathematics, McGill University, Montreal, Quebec, Canada, H3A 2T5 

Abstract 

We consider a model for thin liquid films in a rotating cylinder in the small surface tension limit. 
Using dynamical system methods, we show that the continuum of increasing shock solutions persists 
in the small surface tension limit, whereas the continuum of decreasing shock solutions terminates at 
the limit. Using delicate numerical computations, we show that the existence curves of regularized 
shock solutions on the mass-flux diagram exhibit loops. The number of loops increases and their 
locations move to infinity as the surface tension parameter decreases to zero. If n is the number of 
loops in the mass-flux diagram with 2n + l solution branches, we show that n + 1 solution branches 
are stable with respect to small perturbations. 

1 Introduction 

The time evolution of a liquid film spreading over a solid surface under the action of the surface tension 
and viscosity can be described by lubrication models [21 |8l O [15] . These models approximate the full 
Navier-Stokes equations that appear in the study of motion and instabilities of the liquid film dynamics. 
Thin films play an increasingly important role in a wide range of applications, for example, packaging, 
barriers, membranes, sensors, semiconductor devices, and medical implants. Spin coating is one of the 
methods that is widely used to coat uniform thin films onto solid surfaces in a variety of industrial 
applications such as manufacturing of the magnetic and optical discs. In this technology, a liquid drop 
spreads radially due to centrifugal effects from spinning and eventually yields a thin film of uniform 
thickness formed on the solid surface. In experiments, a variety of different types of bifurcations and 
steady states have been observed [71 \T0\ [IB]. 

In this paper we consider the dynamics of a viscous incompressible thin fluid film on the outer surface 
of a horizontal circular cylinder that is rotating around its axis in the presence of a gravitational field. 
The coating flow is generated by viscous forces due to the cylinder's surface motion relative to the 
fluid. There is no temperature gradient, hence the interface does not experience a shear stress. If the 
cylinder is fully coated there is only one free boundary: where the liquid meets the surrounding air. 
Otherwise, there is also a free boundary (or contact line) where the air and liquid meet the cylinder's 
surface. The motion of the liquid film is governed by four physical effects: viscosity, gravity, surface 
tension, and centrifugal forces. These are reflected in the following parameters: 

• R is the radius of the cylinder; 
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• w is the constant rate of its rotation; 

• g is the acceleration due to gravity; 

• 1/ is the kinematic viscosity; 

• p is the fluid density; 

• 0" is the surface tension. 

These parameters yield three independent dimensionless numbers: 



R^uj _g_ pR^io^ 



Re = ^, 7=,^, We = ^, (1.1) 



where Re is the Reynolds number and We is the Weber number. 

Taking the ratio e = h/R as a small parameter, where h is the average thickness of the liquid, we 
consider the limit e ^ such that 

X = ^e3 and fi = -fRee^ (1.2) 

remain finite and nonzero. 

One can model the fiow on a rotating cylinder using the full Navier-Stokes equations for the velocity 
vector u{r,9, z,t), where r € [0,h{6,t)] is the axial variable, 9 G [— 7r,7r] is the angular variable, z G M 
is the variable in the direction of the cylinder, and h{9,t) is the thickness of the fluid on the surface 
of the cylinder at time t. The simplifying model that takes into account the small surface tension and 
the gravitational force was considered in a number of works including Pukhnachev [13] and O'Brien 
[5]. This model is written in the form, 

dth + de[h- \ph^ cos{9)\ + \xde [h^^ {deh + dlh)] = 0, (1.3) 

where pL and x ^'^^ given in (II. 2p and h{9 + 27r,t) = h{9,t). The model assumes no-slip boundary 
conditions at the liquid/solid interface. A solution to equation (jl.3p is physically relevant if either h 
is strictly positive (the cylinder is fully coated) or h is nonnegative (the cylinder is wet in some region 
and dry in others). 

In a similar context of the syrup rings on a rotating roller, Moffatt [12] neglected the effects of the 
surface tension (i.e., We~^ = = x) and obtained the reduced equation, 

dth + de[h- \ph^ cos(0)] = 0. (1.4) 

The stationary solutions of ()1.3p are given by the 27r-periodic solutions of the third-order differential 
equation, 

h - \ph^ cos{9) + \xh^ {deh + d^h) = Q, (1.5) 

where Q is the constant that corresponds physically to flux of the liquid through the fllm cross section. 

Moffatt [12J and O'Brien & Gath [6j considered stationary solutions of (|1.5p with x = 0- Besides 
smooth periodic solutions for small values of Q, there are two continua of shock solutions for a critical 
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value of (5 = Qc; one is associated with the increasing shocks and the other one is associated with the 
decreasing shocks. Both famihes of solutions are parameterized by the integral 



which has the physical meaning of the mass of the liquid. The solution branches can be plotted on the 
parameter plane (M, Q) for fixed values of /i and X) which we term as the mass-flux diagram. 

Physical arguments convinced the authors of [6] that the decreasing shocks are unstable and cannot 
be observed in the rotational cylinder if the surface tension effects are included with x 7^ 0- The 
singular perturbation theory of small x was recently considered by Benilov et al. [3] , where the authors 
showed that the decreasing shocks do not exist for small positive x- ttis other hand, increasing 
shocks become regularized for small x > ^-iid asymptotic arguments complemented by the numerical 
approximations were developed in [3] to predict spectral stability of stationary solutions with regularized 
increasing shocks. 



Pukhnachev [14] proved the existence and uniqueness of the steady states in the differential equation 

()1.5p if X and Q are not too large. Karabut constructed two branches of steady states in the 
opposite limits of large x- Numerical approximations in Benilov et al. [4J (Figure 14) showed that the 
mass-flux diagram may become more complicated for small x and large M and may include a loop 
near the value Q = Qc- Three solutions coexist for a fixed M if the loop is present. Stability of these 
solutions was not studied in P]. 

It is the purpose of this work to continue, improve, and clarify the preliminary results of Benilov 
et al. OS]. In particular, we develop the dynamical system methods to prove that the family of 
increasing shocks persists with respect to x 7^ 0, whereas the family of decreasing shocks terminate at 
X = 0. We develop a delicate numerical approximation of the steady solutions of equation (jl.Sp to show 
that the number of loops on the mass-flux diagram increases when x is reduced to zero and the location 
of these loops go to infinity. We also apply numerical approximations of eigenvalues of the linearized 
time evolution associated with the lubrication model ()1.3p and show that if n is the number of loops 
on the mass-flux diagram, then n + 1 solution branches are stable with respect to small perturbations. 

The article is organized as follows. Section 2 presents results of the geometric theory on persistence 
of increasing regularized shocks. Section 3 presents numerical results on multi- valued loops in the 
mass-flux diagram. Section 4 provides a summary and discusses open questions. 

Acknowledgments: This work was carried out during the Fields-MITACS Undergraduate Summer 
Research Program in 2010. The authors thank A. Kulyk for collaborations during this program. 

2 Geometric theory of regularized shocks 

We shall study here asymptotic solutions of the steady-state equation (jl.Sp . Using the transformation 
h = Q = Qfi'^^"^, and x = e^^^^ and dropping the tilde sign for h and Q, we obtain the 

third-order differential equation 







3{h - Q) 



e e (-7r,7r), 



(2.1) 



where e > and Q > are parameters of the problem. 
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Solutions of the limiting problem 

FQ{h):=^^^^^ = cos{9), 9e{-7T,7r), (2.2) 

depends on the value of the flux Q > 0. If Q G (O, |), the 27r-periodic solution h{9) is unique [12] . 

Let us denote the smallest roots of FQ{h) = ±1 for Q € (O, |) by h± such that /i_ < hj^. The unique 
solution of the limiting problem (j2.2p satisfies 



h{-9) = h{9): h'{9)>0, 9e{-TT,0), (2.3) 

with h{±7r) = h{0) = h+, and h'{±7r) = h'{0) = 0. 

For Q > |, no solution h{9) exists because max FQ{h) < 1. 



For Q = Q* = |, there is a unique continuous solution h{9) with properties (12. 3p and h{±TT) = h-, 
h(0) = h+ = l, h'i-n) = 0, and lim h'(9) = \. 

Besides this continuous solution at Q = = |, there exists two symmetric families of shock 
solutions with a jump discontinuity at either 9 = 9q ov 9 = —9o, where ^ (O, §) is a continuous 
parameter. 

Let us denote the two simple zeros of Fq^iji) = cos{9q) € (0,1) by H± such that < H^. The 
increasing shock is centered at 9 = —9o and satisfies 

h'{9)>0, 0G(-vr,-0o), h'{9)<0, 9 e {-9o,7t), (2.4) 

with h{±TT) = lim h{9) = lim h{9) = H+, and /i'(±7r) = 0. 

0^00-0 9^00+0 



Using the symmetry the limiting problem (j2.2p with respect to reflection 9^—9, the decreasing 
shock can be constructed using the reflection. It is then centered at 9 = 9q. 

The net mass M deflned by ()1.6p is a one-to-one increasing function of Q for Q € (0,(5*) with 
lim M = and lim M = M* for some M* < oo whereas the two families of shock solutions for 

Q = correspond to the values of M G (M*, oo). 

To consider the persistence of the two shock solutions with respect to parameter e, we shall zoom 
the coordinate 9 near zb^o by the transformation 

h{9) = H{x), x=^-^. (2.5) 

The new function H[x) satisfies a new version of the third-order differential equation 

^+.^/^^ = cos(±^o + eV3.)-F(^), xG^-^,-^j, (2.6) 

where H{x) is a periodic function with period 27re~^/'^ and 

F(//)=Fq.(//) = ^^. 



4 



The limiting problem at e = becomes now the autonomous equation, 



cos{eo) - F{H), xeR. (2.7) 



Recall that H± are zeros of F{H) = cos(0o) with ordering H_ < H^. The increasing shock h{6) 
corresponds to a heteroclinic orbit H(x) of the limiting problem (j2.7p satisfying the boundary conditions 

lim H{x) = H±. (2.8) 

X— >itoo 

The decreasing shock corresponds to a heteroclinic orbit with the boundary conditions 

lim H{x) = H^. (2.9) 

x— >itoo 

Linearization of the limiting equation (|2.7|) near the equilibrium states H± gives 

d^H+ 



dx^ 



-F'{H±)H±, (2.10) 



where F'(H_) > and F'{H^) < 0. Therefore, the equilibrium state H_ has a two-dimensional 
unstable manifold W'^{H-) and a one- dimensional stable manifold W'^{H^), whereas the equilibrium 
state has a one- dimensional unstable manifold VF"(i?+) and a two-dimensional stable manifold 
W'{H+). 

Intersection of the two-dimensional manifolds in W'^{H^) n W^{H^) is transverse in the space M^. 
Hence, a homoclinic orbit satisfying the boundary conditions (|2.8p exists generally and persists under 
the perturbation. 

Intersection of the one-dimensional manifolds in W'^{H^) fi W'^{H-.) is non-transverse in the space 
M^. Hence, a homoclinic orbit satisfying the boundary conditions (j2.9p does not exist generally and 
does not persist under the perturbation. Moreover, it was shown in [3j using sign-definite integral 
quantities that no solution of the limiting equation (|2.7p with the boundary conditions (j2.9p exists. 



As a result, the geometric theory implies that the increasing shock ()2.3p centered at 9 = —6q persists 
as a smooth solution h{9) of the third-order equation (12. ip for any small e > 0, whereas the decreasing 
shock centered at ^ = does not persist in the third-order equation ()2.ip for any small e > 0. 

This conclusion holds for any fixed M > M^,. It does not exclude, however, a possibility of a 
complicated branching behavior in the solutions of the third-order equation ()2.ip that can come from 
M = oo for small values of e > 0. We shall consider the construction of solutions of the third-order 
equation (12. ip numerically. 



3 Numerical analysis 

We shall construct numerical approximations of solutions of the third-order diff'erential equation (12. ip . 
The numerical approximations were generated using a custom-written turning-point algorithm and im- 
plemented in MATLAB. Solutions were found using Newton- Rap hson iterations using Fourier spectral 
differentiation matrices with 256, 512, and 1024 Fourier modes. The mass-fiux diagram was generated 
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with parameter continuation of Q or M, as decided by the algorithm. A convergent solution was defined 
numerically if the (n + 1)*^ iteration hn+i{0) satisfied 

sup |/i„+i(0)-/i„(0)| <5x 10-^ 

0G[-7r,7r] 

Eventually non-convergent solutions were found and the parameter continuation failed at this 'turning 
point'. To resolve the mass-flux diagram near the turning point, the following algorithm was imple- 
mented: first, a reference point was identified in k steps behind the turning point (for our simulations, 
we generally chose k = h). Next, the convergence of the points making up a half-circle centered at the 
turning point was checked. The orientation of the half-circle was chosen to be facing away from the 
direction of the current parametrization (i.e. if we are increasing along the vertical axis then the lower 
half of the circle would be chosen). A vector was then drawn from the turning point to the convergent 
point with the largest distance from the reference point. Finally, the direction of new parametriza- 
tion was chosen from the largest component of this vector. We found that this algorithm successfully 
navigated the loops in the mass-flux diagram. 

Figure [1] shows the mass-flux diagram of stationary solutions for four values of e. For e = 0.005 
(dashed curve), we see no loops in the mass-flux diagram. For each fixed value of mass M, there is 
exactly one value of the flux Q for the stationary solution of (j2.ip . The first loop is formed for e = 0.001 
(light gray). In an interval of values of M, three stationary solutions coexist for three different values 
of Q. Note that this loop was discovered by Benilov et al. [1] (Figure 14). 

Reducing e further, we observe a formation and persistence of the second loop in the mass-flux 
diagram for e = 0.0005 (dark gray) and e = 0.0001 (solid black). Five solutions for different values of Q 
coexist in an interval of values of M. The number of loops keeps increasing as e decreases to zero and 
their location is drifted to large values of M. As e — )• 0, the mass- flux diagram represent an increasing 
curve for M G (0, M^) and a constant level Q = for M > M^, , where = | and ^ 4AA6. The 
limiting picture corresponds to the mass-flux diagram of the solutions of the limiting equation (j2.2p 
that include the smooth solutions for Q < and the shock solutions for Q = Q^,. 

We focus now on two particular examples of the mass-flux diagram with a single loop and a double 
loop. Figure [2] (top) shows the mass-flux diagram for e = 0.001 with a single loop. We can identify 
three solution branches (labeled as S1,S2, and S3) connected at two bifurcation points (labeled as BFl 
and BF2). The other point of intersections of solution branches SI and S2 is not a bifurcation point 
because the two solutions for the same value of M and Q remain distinguishable into two different 
solutions. For M = 10.5, we compute the solution profiles and show them on Figure [2] (bottom). 
Although similar in their shapes, the three steady state solutions are clearly distinct. The peaks of the 
solutions are located for > vr, or equivalently for < 0, thanks to the 27r-periodicity of the solutions. 
They correspond to the increasing shock solution as e ^ located at = —9q < 0. Oscillations, which 
are visible on both sides of the shock are attributed to complex eigenvalues of the linearized equation 
(j2.10p after the scaling transformation (j2.5p and the limit e ^ 0. 

Because of multiple steady-state solutions with the same physical parameter of the mass M, we 
anticipate that they may have different stability properties. Therefore, we examine eigenvalues of the 
linearized equation 



f-h^ cos(0)/ + eh\dgh + dlh)f + ^eh^idef + dif) 



0, (3.1) 



where h{6) is a 27r-periodic steady-state solution and f(9) is a 27r-periodic perturbation to the steady 
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Figure 1: The mass (M) versus the flux (Q) of the steady state solutions of a thin film on the interior 
of a rotating cylinder for various values of e: 0.005 (dashed), 0.001 (light gray), 0.0005 (dark gray), 
and 0.0001 (black). 

state with the growth rate A € C. If Re(A) > 0, the perturbation /(^)e^* grows on the background of 
h{9) and induces instability of the steady-state solution. Expressing h{9) from the third-order equation 
(j2.ip . we can rewrite the spectral problem (13. ip in the equivalent form, 



Note that there is always a zero eigenvalue in the spectral problem (j3.2p . Indeed, /o = ^ is an 
eigenfunction for A = because Q is a free parameter of the solution h. We will now show that the 
zero eigenvalue is simple. First, since the Jacobian operator in the iteration algorithm was found to be 
invertible, the operator Lq is invertible, where 



Since L = OoLq and 8$ has a one-dimensional kernel, the operator L has at most one eigenvector in 
the kernel, which is /o = 

Next, we consider the adjoint spectral problem, 




(3.2) 





(3.3) 



It is clear that go = 1 is the adjoint eigenfunction for A = and that 
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Figure 2: (Top) A segment of the mass-flux diagram for e = 0.001. Three branches are indicated 
between two bifurcation points (labeled BFl and BF2). The three branches are shown by solid gray 
line (SI), solid black line (S2), and dashed line (S3). The dashed gray line shows the value of mass 
M = 10.5. (Bottom) Three steady state solutions with M = 10.5 and e = 0.001: SI (solid gray line, 
Q = 0.6601), S2 (solid black line, Q = 0.6652), and S3 (dashed line, Q = 0.6686). 
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Figure 3: The real part of the smallest eigenvalues A of the spectral problem (j3.ip for e = 0.001 and 
M = 10.5. 



is nonzero at all values of Q but the bifurcation points BFl and BF2. By Fredholm' theory for isolated 
eigenvalue, this fact implies that A = is a simple eigenvalue for all values of Q but the bifurcation 
points BFl and BF2. 

We use a numerical method based on building a matrix representation of the differential operator 
L acting on /. To do so, we discretize the space and approximate the derivatives using the Fourier 
spectral method. The eigenvalue problem is then solved using the MATLAB function eig. 

The real part of the smallest eigenvalues A of the spectral problem (13. 2p is shown on Figure [3] for 
e = 0.001 and M = 10.5. All other eigenvalues have larger negative real parts. Between the two 
bifurcation points BFl and BF2, one eigenvalue crosses zero and becomes unstable, revealing a saddle- 
node bifurcation at points BFl and BF2. Branch S3 (between BFl and BF2) is unstable with exactly 
one real positive eigenvalue A. Other two branches SI and S2 are stable with all but one zero eigenvalue 
having negative real parts. It follows from Figure [2] (bottom) that the profile of the steady-state solution 
at branch S3 is squeezed between those of the solutions at branches SI and S2. Small perturbations 
of the middle tail of the unstable solution is expected to grow towards the upper or lower tails of the 
other two stable solutions. Particular values of A for branches SI, S2, and S3 are given in Table 1. 



solution branch 


real A 


complex A 


SI, Q = 0.6601 


-0.32 


-0.44±i0.79, -0.72±il.50 


S2, Q = 0.6652 


-1.34 


-0.75±i0.44, -1.31±i0.97 


S3, Q = 0.6686 


0.28 


-0.80zbi0.ll, -1.18±il.07 



Table 1: Smallest nonzero eigenvalues of the spectral problem (j3.ip for e = 0.001 and M = 10.5. 



Figure U] (top) shows the mass-flux diagram for e = 0.00039 with two loops. We can identify five 
solution branches (labeled as S1,S2,S3,S4, and S5) connected at four bifurcation points (labeled as 
BFl, BF2, BF3, and BF4). For M = 12.93, we compute the solution profiles and show them on 
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Figure 4: (Top) A segment of the mass-flux diagram for e = 0.00039. Five branches are indicated 
between four bifurcation points (labeled BFl, BF2, BF3, and BF4 and shown by circles). The five 
branches are shown by solid light gray line (SI), solid dark gray line (S2), dashed line (S3), dashed gray 
line (S4), and solid black line (S5). The dashed gray line shows the value of mass M = 12.93. (Bottom) 
Five steady state solutions with M = 12.93 and e = 0.00039: SI (sohd light gray Hue, Q = 0.6638), S2 
(solid dark gray line, Q = 0.6648), S3 (dashed line, Q = 0.6653), S4 (dashed gray line, Q = 0.6661), 
S5 (solid black line, Q = 0.6666). 
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Figure 5: The real part of the smallest eigenvalues A of the spectral problem (j3.ip for e = 0.00039 and 
M = 12.93. 

Figure m (bottom). Although similar in their shapes, the five steady state solutions are clearly distinct. 
Properties of these solutions resemble those on Figure [2j In particular, multiple steady states exist for 
a fixed mass M and e and can be identified by their flux values Q. The five steady-state solutions are 
almost identical with the most visible deviation in their tails and peak heights. 

Comparison between Figures [2] and H] shows that how as e decreases, more loops on the mass-flux 
diagram form and the number of steady states increases. 

Although the five steady-state solutions exist mathematically, whether or not they could exist phys- 
ically depends on their stability. The real part of the smallest eigenvalues A of the spectral problem 
()3.ip is shown on Figure [5] for e = 0.00039 and M = 12.93. All other eigenvalues have larger negative 
real parts. Between the bifurcation points BFl and BF2, BF2 and BF3, BF3 and BF4, one eigenvalue 
crosses zero, revealing four saddle-node bifurcation at these points. Branch S3 (between BFl and 
BF2) is unstable with exactly one real positive eigenvalue A. Branch S4 between BF3 and BF4 is also 
unstable with exactly one real positive eigenvalue. Other three branches SI, S2, and S5 are stable with 
all but one zero eigenvalue having negative real parts. Again, we point readers to Figure U] (bottom) 
that shows how tails of unstable solutions S3 and S4 are aliased by the tails of stable solutions SI, S2, 
and S5. Particular values of A for branches SI, S2, S3, S4, and S5 are given in Table 2. 



branch 


real A 


complex A 


SI, Q = 0.6638 


-0.28, -3.76 


-0.42 ± i0.70, -0.65 ± il.41, -0.90 ± i2.25 


S2, Q = 0.6648 


-0.07, -0.65 


-1.06±i0.82, -1.33 ±fl.82, -1.62 ±i2.86 


S3, Q = 0.6653 


0.30, -0.81 


-1.17±i0.55, -1.50 ±il.56, -1.84±i2.68 


S4, Q = 0.6661 


0.05, -0.69 


-1.03 ± iO.96, -1.25 ± il.94, -1.52 ± i2M 


S5, Q = 0.6666 


-0.69, -2.57 


-1.33±i0.46, -1.54±il.46, -2.31±il.32 



Table 2: Smallest nonzero eigenvalues of the spectral problem (j3.ip for e = 0.00039 and M = 12.93. 
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4 Conclusion 



We have thus explored the behavior of a thin hquid film in a rotating cylinder accounting for surface 
tension and gravity. In particular, we showed that regularized increasing shock solutions persist under 
small surface tension. These shock solutions were then visualized for a wide range of parameter values 
by using numerical discretizations on an uniform grid. We have identified a number of solution branches 
on the mass-flux diagram and have shown numerically that the number of branches increases and the 
location of shocks move to infinity as the surface tension decreases to zero. 

We conclude by listing a number of open questions for further studies. First, it is suggested by the 
numerical computations that the number of solution branches goes to infinity as e — )• but computations 
become difficult and unreachable for e < 10~^. Second, the steady state are expected to persist with 
respect to small inclinations of the cylinder [1] but we do not include inclined cylinders in this work. 
Finally, numerical discretizations on the adaptive (variable) grid can be developed further to resolve 
better the regularized shock solutions near the shock location. 
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